####
###roc绘制
setwd("E:\\1.Ovarian_Treg\\1.picture\\7.ROC")
#引用包
library(survival)
library(survminer)
library(timeROC)

riskFile="nomoRisk.txt"     #风险文件
cliFile="clinical.txt"      #临床数据文件

#读取风险输入文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)
risk=risk[,c("futime", "fustat", "riskScore","Nomogram")]
cli=read.table(cliFile, header=T, sep="\t", check.names=F)
cli=cli[!duplicated(cli$id),]
row.names(cli)=cli$id
cli=cli[,-1]
#合并数据
cli=cli[,-c(1,2)]
samSample=intersect(row.names(risk), row.names(cli))
risk1=risk[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(risk1, cli)

#定义颜色
bioCol=ggsci::pal_npg("nrc")(10)[-2]
#rt$N=as.numeric(rt$N)
#rt=rt[,-10]
######绘制临床的ROC曲线######
predictTime=1     #定义预测年限
aucText=c()
pdf(file="ALLcliROC-1.pdf", width=5, height=5)
#绘制风险得分的ROC曲线
i=3
ROC_rt=timeROC(T=risk$futime,
               delta=risk$fustat,
               marker=risk$riskScore, cause=1,
               weighting='aalen',
               times=c(predictTime),ROC=TRUE)
plot(ROC_rt, time=predictTime, col=bioCol[i-2], title=FALSE, lwd=5)
aucText=c(paste0("Risk", ", AUC=", sprintf("%.3f",ROC_rt$AUC[2])))
abline(0,1)
#rt=rt[,-4]
#unique(rt$Gender)
#对临床数据进行循环，绘制临床数据的ROC曲线
for(i in 4:ncol(rt)){
  ROC_rt=timeROC(T=rt$futime,
                 delta=rt$fustat,
                 marker=rt[,i], cause=1,
                 weighting='aalen',
                 times=c(predictTime),ROC=TRUE)
  plot(ROC_rt, time=predictTime, col=bioCol[i-2], title=FALSE, lwd=5, add=TRUE)
  aucText=c(aucText, paste0(colnames(rt)[i],", AUC=",sprintf("%.3f",ROC_rt$AUC[2])))
}
#绘制图例，得到ROC曲线下的面积
legend("bottomright", aucText,lwd=2,bty="n",col=bioCol[1:(ncol(rt)-1)])
dev.off()

setwd("E:\\1.Ovarian_Treg\\1.picture\\7.ROC")
#引用包
library(dplyr)
library(survival)
library(rms)
library(pec)
riskFile="nomoRisk.txt"     #风险文件
cliFile="clinical.txt"      #临床数据文件

#读取风险文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)
risk=risk[,c("futime", "fustat", "riskScore","Nomogram")]

#读取临床数据文件
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)

#合并数据
samSample=intersect(row.names(risk), row.names(cli))
risk1=risk[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(risk1, cli[,-(1:2)])

#定义颜色
bioCol=ggsci::pal_npg("nrc")(10)

#计算C-index值
riskScore=cph(Surv(futime,fustat)~riskScore, data=rt, surv=TRUE)
Age=cph(Surv(futime,fustat)~Age, data=rt, surv=TRUE)
#Gender=cph(Surv(futime,fustat)~Gender, data=rt, surv=TRUE)

Stage=cph(Surv(futime,fustat)~Stage, data=rt, surv=TRUE)
Nomogram=cph(Surv(futime,fustat)~Nomogram, data=rt, surv=TRUE)
# T=cph(Surv(futime,fustat)~T, data=rt, surv=TRUE)
# M=cph(Surv(futime,fustat)~M, data=rt, surv=TRUE)
# N=cph(Surv(futime,fustat)~N, data=rt, surv=TRUE)
c_index  <- cindex(list("Risk score"=riskScore, 
                        "Nomogram"=Nomogram, 
                        "Age"=Age,
                        #"Gender"=Gender,
                        # "T"=T,
                        # "M"=M,
                        # "N"=N,
                        "Stage"=Stage),
                   formula=Surv(futime,fustat)~ .,
                   data=rt,
                   eval.times=seq(0,6,1),
                   splitMethod="bootcv",
                   B=1000
)
#输出图形
pdf(file="C-index.pdf", width=5, height=4.5)
plot(c_index, 
     xlim=c(3,6), ylim=c(0.25,1.0), 
     col=bioCol, xlab="Time (years)",
     legend.x=5, legend.y=1.1, legend.cex=0.8)
dev.off()

